;;;;;;---------- apply Stream to Jacobi and Gauss-Siedel methods

(load "iterative-methods.scm")
(load "../../../libs-for-chezscheme/stream.ss")
(import (stream))

;;;;;--------No acceleration--------------------
;;---- Jacobi method
(define (Jacobi-stream-from A b n x)
  (stream-cons x
               (Jacobi-stream-from  A b n (nextstate-Jacobi A b x n))))
;;-----Gauss-Siedel method
(define (GaussSiedel-stream-from A b n x)
  (stream-cons x
               (GaussSiedel-stream-from A b n (nextstate-GaussSiedel A b x n))))


;;;;;-----------Aitken acceleration--------------------
(define (Aitken-acceleration strm)
  (let* ([s0 (stream-ref strm 0)]
         [s1 (stream-ref strm 1)]
         [s2 (stream-ref strm 2)]
         [news (vector-map (lambda (e0 e1 e2)
                             (if (zero? (+ e0 e2 (* -2 e1)));;denominate is 0
                                 e2
                                 (- e2 (/ (expt (- e2 e1) 2)
                                          (+ e0 e2 (* -2 e1))))))
                           s0 s1 s2)])
    (stream-cons news (Aitken-acceleration (stream-cdr strm)))))
(define (Aitken-Jacobi-strm A b n x0)
  (Aitken-acceleration (Jacobi-stream-from A b n x0)))
(define (Aitken-GaussSiedel-strm A b n x0)
  (Aitken-acceleration (GaussSiedel-stream-from A b n x0)))


;;;;;-------------super acceleration----------------------
(define (stream-of-streams trans s)
  (stream-cons s (stream-of-streams trans (trans s))))
(define (super-accleration acceler s)
  (stream-map stream-car (stream-of-streams acceler s)))
;;^^^^ above method is general

;;note: accleration should wait some steps
(define (super-Jacobi-strm A b n x0)
  (super-accleration Aitken-acceleration (stream-drop 3 (Jacobi-stream-from A b n x0))))
(define (super-GaussSiedel-strm A b n x0)
  (super-accleration Aitken-acceleration (stream-drop 3 (GaussSiedel-stream-from A b n x0))))



;;------- test
(define (diff sol A b)
  (vector-norm (vector-map - (matrix-*-vector A sol) b)))
(let* ([A (vector '#(0.8  1.   0.9)
                  '#(0.7 0.5   1)
                  '#(0.6 0.5   1))];;这个矩阵可谓是最糟糕的了, 但GaussSiedel+Aitken仍成功
       [b (vector 1. 2. 3.)]
       [n 3]
       [x0 (vector 0. 0. 0.)])
  (newline)
  (display "jacobi orgin: \n")
  (for-each (lambda (i)
              (display (diff (stream-ref (Jacobi-stream-from A b n x0) i) A b))
              (newline))
            (range 0 15 5))
  (display "gs orgin: \n")
  (for-each (lambda (i)
              (display (diff (stream-ref (GaussSiedel-stream-from A b n x0) i) A b))
              (newline))
            (range 0 15 5))
  (display "aikten-jacobi: \n")
  (for-each (lambda (i)
              (display (diff (stream-ref (Aitken-Jacobi-strm A b n x0) i) A b))
              (newline))
            (range 0 15 5))
  (display "aitken-gs: \n")
  (for-each (lambda (i)
              (display (diff (stream-ref (Aitken-GaussSiedel-strm A b n x0) i) A b))
              (newline))
            (range 0 15 5))
  (display "super-jacobi: \n")
  (for-each (lambda (i)
              (display (diff (stream-ref (super-Jacobi-strm A b n x0) i) A b))
              (newline))
            (range 0 15 5))
  (display "super-gs: \n")
  (for-each (lambda (i)
              (display (diff (stream-ref (super-GaussSiedel-strm A b n x0) i) A b))
              (newline))
            (range 0 15 5))
  )

